Problème torsion d’arbres#

Librairies#

Hide code cell content
# Importation des librairies
import dolfinx  # Module principal de FEniCSx pour le calcul par éléments finis
from dolfinx import mesh, fem, plot, io, default_scalar_type  # Sous-modules FEniCSx essentiels
from dolfinx.fem.petsc import LinearProblem  # Résolution de systèmes linéaires
from mpi4py import MPI  # Interface pour le calcul parallèle
import ufl  # Unified Form Language pour les formulations variationnelles
import numpy as np  # Calcul numérique efficace
import matplotlib.pyplot as plt
import math  # Module de fonctions mathématiques standard de Python
import pyvista as pv  # Visualisation 3D scientifique
import gmsh  # Générateur de maillage 3D
import meshio  # Lecture/écriture de différents formats de maillage
from dolfinx.io import XDMFFile  # Gestion des fichiers XDMF pour données volumineuses
import dolfinx.fem as fem  # Module FEniCSx pour la méthode des éléments finis
import dolfinx.mesh as mesh  # Module FEniCSx pour la gestion avancée des maillages
from petsc4py import PETSc  # Suite de solveurs numériques pour problèmes à grande échelle
import panel as pn
from myst_nb import glue

Géométrie#

Définition de la géométrie et du maillage#

Hide code cell content
# Initialisation de GMSH
gmsh.initialize()
gmsh.model.add("ellipse")

# Définition des paramètres géométriques et de maillage
a = 0.2  # Demi-grand axe de l'ellipse en mètres
b = 0.2  # Demi-petit axe de l'ellipse en mètres
h = 1.0  # Hauteur du cylindre elliptique en mètres
# Calcul du rayon équivalent du cylindre elliptique
R = 1/2 * np.sqrt(a**2 + b**2)

# Création de la géométrie : ellipse de base
ellipse = gmsh.model.occ.addEllipse(0, 0, 0, a, b)

# Création d'une boucle fermée à partir de l'ellipse
curve_loop = gmsh.model.occ.addCurveLoop([ellipse])

# Création d'une surface plane à partir de la boucle
surface = gmsh.model.occ.addPlaneSurface([curve_loop])

# Extrusion de la surface pour créer un cylindre elliptique
gmsh.model.occ.extrude([(2, surface)], 0, 0, h)

# Synchronisation du modèle géométrique
gmsh.model.occ.synchronize()

# Configuration des paramètres de maillage
lc = 0.02  # Taille caractéristique des éléments du maillage en mètres
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)

# Génération du maillage volumique (3D)
gmsh.model.mesh.generate(3)

# Sauvegarde du maillage au format .msh
gmsh.write("ellipse.msh")

# Finalisation de GMSH
gmsh.finalize()

# Lecture du fichier .msh avec meshio
msh = meshio.read("ellipse.msh")

# Conversion et sauvegarde du maillage au format XDMF (compatible avec FEniCSx)
meshio.write("ellipse.xdmf", meshio.Mesh(points=msh.points, cells={"tetra": msh.cells_dict.get("tetra", [])}))

# Lecture du maillage XDMF avec FEniCSx
with XDMFFile(MPI.COMM_WORLD, "ellipse.xdmf", "r") as xdmf_file:
    domain = xdmf_file.read_mesh(name="Grid")
    domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)

# Affichage des informations sur le maillage
print("Maillage créé et importé avec succès.")
print(f"Nombre de sommets : {domain.topology.index_map(domain.topology.dim).size_global}")
print(f"Nombre d'éléments : {domain.topology.index_map(domain.topology.dim-1).size_global}")

Visualisation de la géométrie#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=700, height=400)
interactive_panel
from myst_nb import glue
glue("img", interactive_panel,display=False)

Configuration du problème#

import ipywidgets as widgets
from IPython.display import display

# Exemple de slider pour changer le rayon de la poutre
rayon = widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01, description="Rayon")
display(rayon)

# Code Fenics ou autre qui dépend de "rayon.value"
import ipywidgets as widgets
from IPython.display import display

# Exemple de slider pour changer le rayon de la poutre
rayon = widgets.FloatSlider(value=0.05, min=0.01, max=0.1, step=0.01, description="Rayon")
display(rayon)

# Code Fenics ou autre qui dépend de "rayon.value"

Définition de l’espace fonctionnel#

# Exemple de dropdown pour choisir un matériau
materiau = widgets.Dropdown(
    options=['Acier', 'Aluminium', 'Plastique'],
    value='Acier',
    description='Matériau:',
)

# Fonction pour mettre à jour le slider en fonction du matériau choisi
def update_slider(change):
    if change['new'] == 'Acier':
        rayon.value = 0.05  # Valeur par défaut pour l'acier
    elif change['new'] == 'Aluminium':
        rayon.value = 0.04  # Valeur par défaut pour l'aluminium
    else:
        rayon.value = 0.03  # Valeur par défaut pour le plastique

materiau.observe(update_slider, names='value')
imk=materiau.observe(update_slider, names='value')

display(materiau)
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(imk.ren_win, width=700, height=400)
from myst_nb import glue
glue("img", interactive_panel)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[101], line 3
      1 # Utiliser Panel pour créer un rendu interactif
      2 pn.extension('vtk')
----> 3 interactive_panel = pn.pane.VTK(imk.ren_win, width=700, height=400)
      4 from myst_nb import glue
      5 glue("img", interactive_panel)

AttributeError: 'NoneType' object has no attribute 'ren_win'

Définition des frontière du domaine#

Hide code cell content
# Définition des fonctions pour identifier les différentes faces du cylindre
def down_face(x):
    """Identifie la face inférieure du cylindre (z = 0)"""
    return np.isclose(x[2], 0)

def top_face(x):
    """Identifie la face supérieure du cylindre (z = h)"""
    return np.isclose(x[2], h)

def lateral_face(x):
    """
    Identifie la surface latérale du cylindre elliptique
    Utilise l'équation de l'ellipse : x^2/a^2 + y^2/b^2 = 1
    """
    tolerance = 1e-5  # Tolérance pour la comparaison numérique
    return (np.isclose((x[0]**2 / a**2 + x[1]**2 / b**2), 1.0, atol=tolerance)) & (0 <= x[2]) & (x[2] <= h)

# Dimension des facettes (2D pour un domaine 3D)
fdim = domain.topology.dim - 1

# Localisation des entités (facettes) correspondant à chaque face
Sigma_l = mesh.locate_entities(domain, fdim, lateral_face)
Sigma_0 = mesh.locate_entities(domain, fdim, down_face)
Sigma_h = mesh.locate_entities(domain, fdim, top_face)

# Préparation des marqueurs de facettes pour les conditions aux limites
marked_facets = np.hstack([Sigma_h])  # Utilisation uniquement de la face supérieure
marked_values = np.hstack([np.full_like(Sigma_h, 1)])  # Marqueur 1 pour la face supérieure

# Tri et création des tags de maillage
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, domain.topology.dim - 1, 
                                  marked_facets[sorted_facets], 
                                  marked_values[sorted_facets])

# Définition de la mesure pour les facettes marquées
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)

Visualisation des frontière du domaine#

Hide code cell outputs
# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(pl.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Conditions aux limites#

# Définition des conditions aux limites de Dirichlet (déplacement nul)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)

# Application des conditions aux limites sur chaque face
bc1 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_0), V)  # Face inférieure
bc2 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_h), V)  # Face supérieure
bc3 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, Sigma_l), V)  # Surface latérale

Propriétés matériau#

# Définition des constantes élastiques du matériau
mu = 80e9      # Module de cisaillement (G) en Pa
lambda_ = 120e9  # Premier paramètre de Lamé (λ) en Pa

Définition des chargements#

α = 20.0 m^-1
C = 1005309649.1487346 kg.m^2/s

Résolution#

Hide code cell content
# Définition des tenseurs de déformation et de contrainte
def epsilon(u):
    """Tenseur de déformation linéarisé"""
    return ufl.sym(ufl.grad(u))

def sigma(u):
    """Tenseur des contraintes (loi de Hooke)"""
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

# Définition des fonctions d'essai et de test
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Définition de la force volumique (nulle dans ce cas)
f = fem.Constant(domain, default_scalar_type((0, 0, 0)))

# Formulation variationnelle du problème
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx  # Forme bilinéaire
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds(1)  # Forme linéaire

# Définition et résolution du problème linéaire
problem = LinearProblem(a, L, bcs=[bc1], 
                        petsc_options={"ksp_type": "cg", "pc_type": "jacobi"})
uh = problem.solve()

Visualisation des résultats#

Visualisation des déplacements#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Visualisation des rotations principales#

Visualisation des déplacements amplifiés#

Visualisation des déplacements normalisés et selon les axes principaux#

Visualisation des composantes du tenseur des déformations#

Visualisation des composantes du tenseur des contraintes#

Calcul des contraintes de von Mises#

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

V_von_mises = fem.functionspace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

Visualisation des contraintes de von Mises#

# Utiliser Panel pour créer un rendu interactif
pn.extension('vtk')
interactive_panel = pn.pane.VTK(p.ren_win, width=600, height=400)
from myst_nb import glue
glue("img", interactive_panel)

Analyses des résultats#

Calcul des valeurs maximales#

Hide code cell content
# Extraction de la valeur maximale de la rotation RZ
max_rotation_z = np.max(np.abs(rotation_z.x.array))

# Extraction de la valeur maximale de la norme de rotation
max_rotation_norm = np.max(rotation_norm.x.array)

# Conversion en degrés pour une interprétation plus intuitive
max_rotation_z_degrees = np.degrees(max_rotation_z)

# Calcul de l'angle de torsion théorique
theoretical_twist = ang * h  # ang est l'angle de torsion par unité de longueur, h est la hauteur du cylindre

# Comparaison pour RZ
relative_error_Z = abs(max_rotation_z - theoretical_twist) / theoretical_twist * 100

# Comparaison pour norme R
relative_error_N = (max_rotation_norm - theta_theoretical) / theta_theoretical * 100
Moment d'inertie polaire J : 6.283185e-04 m^4
Angle de torsion théorique θ : 20.000000 radians
Angle de torsion théorique θ : 1145.92 degrés

Comparaison avec la simulation numérique en RZ:
Rotation Z maximale simulée : 19.888688 radians
Rotation Z maximale simulée : 1139.54 degrés
Erreur relative : 0.56%

Comparaison avec la simulation numérique en norme:
Norme de rotation maximale : 19.901232 radians
Norme de rotation maximale : 1140.26 degrés
Erreur relative : -0.49%

Valeurs maximales des composantes de rotation :
RX max : 124.69 degrés
RX max  : 2.176298 radians
RY max : 122.15 degrés
RY max  : 2.131950 radians